Efficient Computation of Hyperspherical Bessel Functions 



Arthur Kosowsky* 

Department of Physics and Astronomy, Rutgers University, 136 Frelinghuysen Road, Piscataway, New Jersey 08854-8019 

(May 13, 1998) 

Fast and accurate computations of the power spectrum of cosmic microwave background fluctu- 
ations are essential for comparing current and upcoming data sets with the large parameter space 
of viable cosmological models. The most efficient numerical algorithm for power spectrum calcula- 
tion, recently implemented by Seljak and Zaldarriaga, involves integrating sources against spherical 
Bessel functions or, in the cases of a non-flat universe, analogous hyperspherical Bessel functions. 
Evaluation of these special functions usually dominates the computation time in non-flat spatial 
geometries. This paper presents a highly accurate and very fast WKB approximation for comput- 
QQ ' ing hyperspherical Bessel functions which will greatly increase the speed of microwave background 

' power spectrum computations in open and closed universes. 
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I. INTRODUCTION 



The power spectra of temperature and polarization fluctuations of the cosmic microwave background contain a 
rich harvest of cosmological information. From upcoming high-resolution maps of the microwave background, we can 
. hope to determine the basic cosmological parameters to high precision, as well as test the nature of the primordial 

i: perturbations and the mechanism for structure formation in the Universe . Extracting this information will require 

J computationally intensive analysis of the parameter space of cosmological models; a basic, general family of inflation- 

^-H like models requires around ten parameters. Monte Carlo explorations of such a large parameter space will require 
power spectra evaluations for millions of models, so it is imperative that the fastest possible code be available for 
calculating power spectra. 

' Traditionally, microwave background power spectra were computed by expanding the angular dependence of the 
radiation field in multipole moments and then evolving the resulting set of coupled differential equations. For power 
. spectra at sub-degree scales, typically several thousand coupled equations need to be evolved; the CPU time for such 
Oh' codes on current workstations is generally measured in hours. Even scaling up to supercomputer speeds, large Monte 
Carlo calculations are not possible with such codes. 
' A major advance was the realization by Seljak and Zaldarriaga that a formal integral solution to the photon 
^ [ evolution equation allows the solution (in flat space) to be expressed as a source term integrated against spherical 
Ci • Bessel functions [||. The source term still needs to be calculated via coupled evolution equations, but since only the 
^ [ lowest moments contribute to the source term, far fewer equations must be evolved. In essence, the full set of coupled 
equations integrates the Bessel equation, whose solution is already known. The spherical Bessel functions depend 
on radius r and on two parameters, the wavenumber k and the index I; however, in flat space, the wavenumber can 
be scaled with the radial variable so ji{kr) is actually a one-parameter family of functions. Microwave background 
computations generally require values of Z up to a few thousand, and it is computationally feasible to precompute and 
cache the required Bessel functions. Thus the flat-space power spectrum can be evaluated with a greatly increased 
speed using the Seljak-Zaldarriaga algorithm. 

For open or closed universes, the analogous "hyperspherical" Bessel functions depend on the radial coordinate 
and the wavenumber separately, so an inordinate number of functions would need to be precomputed and stored. 
Seljak and Zaldarriaga resort to a numerical integration of the differential equation defining the hyperspherical Bessel 
functions, which is far slower than calling the precomputed functions in the fiat case. This paper calculates the WKB 
approximation to the hyperspherical Bessel functions, which, while not as fast as precomputation, offers a dramatic 
increase in computational speed over other methods. Furthermore, the WKB approximation is highly accurate for all 
but the lowest few values of I. 

The following Section gives a brief review of WKB theory. Section III then gives an overview of the relevant 
properties of the hyperspherical Bessel functions; Section IV derives the WKB approximation for all cases, including 
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the usual spherical Bessel functions in the flat space limit. Comparisons with exact functions demonstrate the 
remarkable accuracy of the approximation. The paper concludes with brief remarks discussing implementation and 
application of the approximations. 



II. REVIEW OF WKB THEORY 



The WKB approximation applies to equations of the Schrodinger form 

t^y" = Q{x)y. 



(1) 



Such equations have exponential behavior in regions where Q{x) > (dissipative regions) and oscillatory behavior 
when Q{x) < (dispersive regions). Points with Q{x) = are called turning points, marking the transition between 
the two behaviors. The WKB approximation (away from the turning points) consists of an exponentiated power series 
in e: 



y{x) = exp 



^ OO 

-^5„(x)e" 



(2) 



This paper considers the usual first-order WKB approximation, which retains only the n = and n = 1 terms of the 
expansion. This approximation is often surprisingly accurate for e as large as 1. For a more detailed exposition of 
WKB theory, see Ref. |]. 

Equations possessing a single turning point are the most straightforward application of WKB theory. Without loss 
of generality, translate the dependent variable so that the turning point is at the origin and Q{x) > for a; > 0. In 
the dissipative region a; > 0, which will be called region I, a simple asymptotic expansion gives the familiar "WKB 
formula" for the exponentially decaying solution 



yjix) C [Q(x)]-i/^ 



exp 



VQ(f)dt 



(3) 



where C is an undetermined normalization constant. In region II in the neighborhood of a; = 0, Q{x) can be replaced 
by its asymptotic expansion Q ~ ax, x —> which converts Eq. (n]) to an Airy equation with solution 



yiiix) - 2C\/^(ae)-i/^Ai (a^/^e-^/a^,^ 







(4) 



where Ai(x) is the first type of Airy function and the prefactor has been determined by an asymptotic match with 
Eq. (^. Finally, in region III with a; < 0, an analytic continuation of Eq. ^) and an asymptotic match to yjj gives 



yjuix) ~ 2C [Q{x)]-^/^sm \- f ^/^Qif)dt + j 



(5) 



Remarkably, there exists a single function which uniformly approximates y{x) and reduces to the above asymptotic 
forms in the various limits: 



y{x) ~ 2^C 



35'o(x) 
2e 

S'o(x) 



1/6 



-1/4 



Ai 



35'o(a;) 
2e 



2/3' 



(6) 



VQ(i)dt. 



This formula is the basis for all approximations presented in this paper. 

For the case of a function with two turning points at x — A and x = B, two single-turning-point solutions can be 
asymptotically matched to form the general solution, which leads to the eigenvalue condition 



x/-Qit)dt 



n + - \ TT 



(7) 



with n a non-negative integer. This condition typically is used to approximate the energy eigenvalues of a particle in 
an arbitrary potential well in quantum mechanics. 
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III. PROPERTIES OF HYPERSPHERICAL BESSEL FUNCTIONS 



It is straightforward to apply these expressions to the case of hyperspherical Bessel functions. A review of their 
relevant properties is presented here; a more detailed exposition of these functions on which the following discussion is 
based is given in Ref. Q|. The hyperspherical Bessel functions $f (f) are radial eigenfunctions of the covariant Laplace 
operator in spherical coordinates, for spaces of constant curvature: 

(V2 + fc2)$/3(^)y,,^(0^^) ^0 (8) 



where the index j3 — \/k? + K with K — H'^{Vl{) — 1) the spatial curvature, and V is the covariant 3-derivative 
associated with the spatial part of the Robertson- Walker metric: 



ds^ = dt^ - R^{t) 



l-Kr^ 



(9) 



In this paper, all physical distances will be expressed in units of the curvature scale, giving K — 1 for a closed 
universe, K = for a critical (flat) universe, and K = —1 for an open universe. The eigenfunctions $f , termed 
"hyperspherical" or "ultraspherical" Bessel functions because they reduce to the usual spherical Bessel functions in 
the case of a flat universe, are very useful because they serve as the analog of Fourier modes for cosmological models 
with non-zero curvature. For a critical density universe with no spatial curvature, V is just the usual gradient operator 
and $f (r) = ji {(3r) = ji (kr) . 

It is convenient to change variables to the radial coordinate x defined by dx — dr j^/X — Kr'^; explicitly, 

f sinx, K = 1] 

r{x) ^{x, K = 0; (10) 




Then it is straightforward to demonstrate from Eq. (g) that the hyperspherical Bessel functions satisfy the Schrodinger 
equation 



dx^ 



r(x)2 



(11) 



where uf (x) = ^{x)^f ix)- The hyperspherical Bessel fimctions are normalized so as to match the normalization of 
ji{kr) in the flat-space limit. For the cases K = Q and K = —1, the momentum variable [3 can have any positive 
value, Eq. (^l]) has a single turning point for x > 0, and the WKB approximation can be applied directly. 

The K — 1 case is more complicated: the spatial sections of the spacetime are compact, resulting in a discrete 
eigenvalue spectrum for the eigenfunctions. It is possible to express the closed-universe functions in terms of associated 
Legendre polynomials as 

$f (X {s\nx)-"^Pp[\]i\cosx); (12) 

this form shows that $^ is periodic in x with period 27r. As with spherical Bessel functions, <I>^ is symmetric 
(antisymmetric) around x = for ^ even (odd); thus the fimction is determined by its values on the interval [0, tt]. 
Requiring <I>f to be single-valued gives the condition 

$f(-cosx) = cos[(/3 - /- l)7r]$f(cosx)- (13) 

Two conclusions follow immediately: [3 must be a positive integer, and the functions arc symmetric (antisymmetric) 
around x = 7r/2 if /3 — / — 1 is even (odd). Thus the value of <i>f (x) must be computed only on the interval [0, 7r/2] 
and can be determined for all other values of x by symmetry. It can be demonstrated that (3=1 and j3 = 2 represent 
gauge modes, not physical perturbations [^, so /3 takes integer values of 3 or larger, and /3 > I follows from Eq. (p^). 

The hyperspherical Bessel functions also satisfy several useful recursion relations which are given in Ref. A 
closed-form expression for $f can be obtained as Z -f 1 derivatives of an elementary function; it is practicable to use 
these exact expressions for evaluating the hyperspherical Bessel functions up to Z = 4 or 5. For precise calculation of 
the functions for larger values of /, recursion techniques can be used for the open case, in analogy with the Miller's 
method evaluation of Bessel functions in Ref. |^]. However, for the closed case2_downwards recursion is not always 
available due to the restriction I < (3, so direct integration of Eq. (O) is better 
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IV. WKB APPROXIMATIONS 



To apply the WKB formalism most effectively, first divide botli sides of Eq. (^l|) by /(/ + 1) to obtain Eq. witli 

1 



Qix) 



r2(x) 



a 



with a = (3e and e = 1/ + 1). Turning points xo ^i'g located where 



(14) 



(15) 



note that the K = Q and K ~ 1 cases each have a single turning point, while the K = 1 case has a single turning 
point within the range [0,7r/2]. Thus Eq. (^ can be applied directly to each case. The WKB approximation is an 
asymptotic series in powers of e, which in this case is the inverse of I. For larger I values, the approximation becomes 
progressively better; it is also better away from the region of the turning point where the various asymptotic solutions 
are matched. As demonstrated below, the approximation is remarkably good even for / = 2 and 3. 

The WKB approximation offers a great increase in numerical speed because the required integrals can be performed 
exactly in terms of elementary functions. For the K = case. 



-rr — a 



1/2 



dt 



1/2 



dt = log 




- Vl - a^x^. 



For the K — ~1 case, defining w = a sinhx. 



sinh^ t 



1/2 



dt 



f 1 



V sinh^ t 



1/2 




dt ^ — tan ^ 
2 



And similarly, for the K — 1 case, defining v = asinx. 



sin^ t 



1/2 



dt = tan 



- 1 



sin^ t 



1/2 



dt = tanh 



1/2- 



1/2- 



■ tan 



2 

' log 



2^/(^2 _ 1)(q,2 _„2) 

2u2 -a^ -1 



TT 

2' 



(16a) 
(16b) 

(17a) 
(17b) 

(18a) 
(18b) 



Care must be taken to use the correct branch of the inverse tangent functions. 
In the closed case, the eigenvalue condition Eq. (j^) becomes 



1 1 

2+7 



(19) 



with n an integer, so for a given exact integer value of /?, the corrected WKB eigenvalue can be obtained by the 
replacement 



1 

81 



1 

16P' 



(20) 



which is sufficiently accurate for a first-order WKB approximation. If the eigenvalue /3 is not corrected when evaluating 
the approximate function, the function or its first derivative will be discontinuous at x = 7r/2. 

The other necessary numerical ingredient is the evaluation of the Airy function in Eq. (ffl). If a fast routine is not 
available, a reasonable approximation is to use the leading asymptotic behavior at large arguments and a Taylor series 
around the origin. With a crossover at = 1.6 and a series including x^^ terms, the residual error in Ai(a;) is at 
the 1% level. Fast routines based on Chebyshev polynomial fits or Fade expansions are readily obtainable in both 
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FIG. 1. The WKB approximation for the usual spherical Bessel functions ji{x), for Z = 2 and / = 5. The exact functions are 
solid lines, the approximate functions dashed lines. For I > 5 the exact and approximate functions are indistinguishable on the 
scale of the plot. 
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Fortran and C. Note that since the Airy function is evaluated for every hyperspherical Bessel function, errors in Airy 
function evaluation translate into systematic errors in $f . Using the simple asymptotic approximation to Ai(x), for 
example, leads to a systematic 1% error when integrating smooth functions against <i>^ independent of I; this error is 
reduced to 0.01% or less when an accurate Airy function evaluation is employed ||]. 

Figure 1 displays the exact and WKB-approximated spherical Bessel functions ji{x) by I — 2 and I = 5. The 
accuracy is very good even for / = 2, with an error of 1.5% at the first peak in the I = 2 case and 0.6% at the first 
peak in the / = 5 case. For higher values of I, the actual and approximated functions are indistinguishable on the 
scale of the plot; by / = 20, the first peak is accurate to 0.05%. As mentioned above, the harmonics for the lowest 
few I values can be evaluated exactly in terms of trigonometric functions. 




FIG. 2. The WKB approximation for the open universe hyperspherical Bessel functions $1 and shown as the dashed 
lines, compared with the solid exact functions. 



Figures 2 and 3 compare the exact and WKB-approximated hyperspherical Bessel functions for open and closed 
universes, respectively. The level of accuracy is essentially the same as in the flat case, with approximations for I > 5 
indistinguishable from the exact values on the scales of the plots. One particular set of closed- universe functions 
always has a significant error: the lowest eigenvalue corresponding to a given I with the unperturbed value f3 = I + I. 
As shown in Figure 4, the approximate functions have discontinuous derivatives at x — i'/2 for all I. The reason is 
that the turning point xo is so close to 7r/2 for the lowest eigenvalue that the solution does not have enough room 
between xo and 7r/2 to attain its asymptotic behavior away from the turning point, so an asymptotic match around 
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FIG. 3. The WKB approximation for the closed universe hyperspherical Bessel functions $2 and $5, shown as the dashed 
lines, compared with the solid exact functions. 
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X = T^I'^ is unsuccessful. This behavior can be circumvented with an exact form for the lowest eigenf unctions: 



{21)1 



(Z + 1)(2Z + 1)!! 



1/2 



sm X 



{K = l) 



27r 



21 + 1 



1/4 



I 



{l + l){2l + l) 



1/2 



7 



m 512/2 



sin' X, 



I — > oo; 



(21) 



the asymptotic expression is good to 0.02% at Z = 2. Fortunately, the second and higher eigenvalues are not greatly 
affected by this problem (see Fig. 3, for example). 




FIG. 4. The WKB approximation for the closed universe hyperspherical Bessel functions $2 and $5, shown as the dashed 
lines, compared with the solid exact functions. The discrepancy near x = 7i"/2 is because the matching solution does not have 
enough space between the two turning points to reach its asymptotic value. This artifact persists in the lowest eigenfunction 
for all values of I. 
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V. CONCLUDING REMARKS 



The above calculations demonstrate that the WKB approximation to hyperspherical Bessel functions is highly 
accurate, with increasing accuracy as I becomes larger. Furthermore, the computation of one of the functions at 
a given value of x requires around ten elementary function calls, plus some arithmetic, so the approximation is 
much faster than evaluation based on recursion methods or integration of Eq. The computation time is also 

independent of in marked contrast to other methods, and is completely stable for any values of I and (3. This 
approximation will substantially speed the computation of microwave background power spectra in open or closed 
universes, potentially by an order of magnitude based on rough timings of the CMBFAST code Q . 

A further speedup in evaluating these functions can be accomplished by precomputing and caching the integrals in 
Eqs. (^6|) to (]l^); each is a function of the two variables x and a. The integrals are very uneventful functions of these 
variables, and accurate interpolation is possible based on a small number of computed values. With this scheme, the 
computational work for a hyperspherical Bessel function is reduced to evaluating a power and a sine or exponential 
in Eq. (||) or (H), an inverse sine or hyperbolic sine to evaluate xo, and a small amount of arithmetic. 

The uniform WKB approximation presented here is useful for evaluating any functions defined by a second-order 
differential equation in Schrodinger form, especially when speed of evaluation is more important than accuracy to 
many significant figures. This is often the case in physics problems, and particularly in the case where a special 
function is integrated against another function. If it is impracticable to cache all needed values of the special function, 
WKB offers a fast, simple, and remarkably accurate alternative. 
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